package scu.maqiang.meshgeneration;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.SRMatrix;

public class CircleMeshBySolveLaplace {
    public static void main(String[] args) {
        Mesh2T3 mesh = new Mesh2T3().square2D(32, 32, T3Type.Left);
        FES2T31 fs = new FES2T31(mesh);
        int dof = fs.GetNdof();
        SRMatrix A = new SRMatrix(dof);
        double[] RHS = new double[dof];
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        ScalarFunc Xbd1 = (xy, label, param) ->  Math.sin(0.5 * Math.PI * xy[0]);
        ScalarFunc Xbd2 = (xy, label, param) ->  Math.cos(0.5 * Math.PI * xy[1]);
        ScalarFunc Xbd3 = (xy, label, param) -> -Math.cos(0.5 * Math.PI * xy[0]);
        ScalarFunc Xbd4 = (xy, label, param) -> -Math.sin(0.5 * Math.PI * xy[1]);
        fs.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
        ScalarFunc[] bdFunc = {Xbd1};
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 1);
        bdFunc[0] = Xbd2;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 2);
        bdFunc[0] = Xbd3;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 3);
        bdFunc[0] = Xbd4;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 4);
        DirectSSolver dss = new DirectSSolver(A);
        dss.Factorize_Cholesky(1);
        double[] CircleX = dss.Solve_Cholesky(RHS, 1);

        ScalarFunc Ybd1 = (xy, label, param) -> -Math.cos(0.5 * Math.PI * xy[0]);
        ScalarFunc Ybd2 = (xy, label, param) ->  Math.sin(0.5 * Math.PI * xy[1]);
        ScalarFunc Ybd3 = (xy, label, param) ->  Math.sin(0.5 * Math.PI * xy[0]);
        ScalarFunc Ybd4 = (xy, label, param) -> -Math.cos(0.5 * Math.PI * xy[1]);
        bdFunc[0] = Ybd1;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 1);
        bdFunc[0] = Ybd2;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 2);
        bdFunc[0] = Ybd3;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 3);
        bdFunc[0] = Ybd4;
        fs.applyBC_MBN(RHS, Direct.All, bdFunc, null, 4);
        double[] CircleY = dss.Solve_Cholesky(RHS, 1);

        mesh.setNodeCoordinates(Direct.X, CircleX);
        mesh.setNodeCoordinates(Direct.Y, CircleY);
        mesh.toTecplot("CircleMesh.dat");

    }
}
